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Abstract 

This  study  explores  the  problem  of  dynamic  stall,  i.e. 
the  stall  of  an  airfoil  undergoing  pitching  motion.  The 
general  equations  of  continuity  and  momentum  are  developed 
for  a  non-inertial  and  unsteady  control  volume.  They  are 
written  in  momentum-integral  form  and  the  boundary  layer  on 
the  pitching  airfoil  is  computed  using  a  modified  von  Karman- 
Pohlhausen  method. 

The  boundary  layer  edge  velocity,  velocity  gradient  and 
time  rate  of  change  of  velocities  required  for  the  step  by 
step  integration  of  the  von  Karman-Pohlhausen  working 
equations  are  obtained  from  the  inviscid  solution.  The 
inviscid  velocity  profile  along  the  surface  of  the  airfoil  is 
obtained  by  conformal  mapping  from  the  velocity  profile 
around  a  rotating  circular  cylinder.  Complex  potential  flow 
theory  is  used  to  obtain  the  velocity  around  the  cylinder. 
The  Kutta  condition  is  continuously  maintained  at  the  point 
mapping  to  the  trailing  edge  of  the  airfoil  for  each  time 
step.  This  way ,  the  flow  is  considered  steady  at  each  time 
step,  but  varies  from  one  time  step  to  the  next  when  the 
angle  of  attack  is  increased.  The  increase  in  stall 
angle  of  attack  is  analyzed  as  a  function  of  a  non- 
dimensional  pitch  rate  (0.5ca/Uw).  Although  the  solution  is 
obtained  primarily  for  a  symmetric  Joukowski  airfoil  of 
thickness  ratio  0.15  (J015),  the  analysis  also  includes 
variations  of  camber,  thickness  and  location  of  the  point  of 
rotation . 


x 


for  the  general  unsteady,  non-inertlal  control  volume  fixed 
on  the  surface  of  our  pitching  airfoil. 

Figure  4  is  an  expanded  view  of  this  control  volume 
showing  all  the  x-components  of  the  forces  acting  on  and 
within  it.  The  control  volume  is  of  unit  depth, 
infinitesimally  narrow  (dx)  and  of  height  h  such  that  its 
upper  side  is  al\.ays  above  the  edge  of  the  boundary  layer. 

Applying  (20)  to  the  c.v.  of  Figure  3,  we  get  the  x- 
component  of  the  momentum  equation  for  our  c.v.: 

-  3  Pdydx  -  xwdx  -  R(asin<j>  +  a2cos$)  +  2av  +  ay 

3xo  o 

-  a2x}  pdydx  =  3_  /h ppu  dxdy  +  lLmTnp  +  3  /^pu  dxdy 

3x  o  1U  3t  o  [  ’ 

Now,  applying  the  continuity  equation  (1)  to  the  same  control 
volume  and  assuming  incompressible  flow  (p  =  constant),  we 
get  from  Figure  5: 

m-j-Qp  =  -3_  /hpu  dydx  (22) 

3x  o 

Using  (22)  in  (21)  yields: 

l_  i_ 

■3 _  /  Pdydx  -  xwdx  /  {  R(asirK)>  +  azcoscJ>)  +  2av  +  ay  -a2}  p  dydx 
Ox  o  o 

=  3_  /hpuudydx  -  U  _3_  /hpudydx  +  3__  /hpu  dydx 

3x  o  3x  o  3t  o 


(23) 


We  assume  that  our  airfoil  is  perfectly  rigid  so  that  the 
quantities  R,  0  ,  <p ,  ip  are  all  constant  while  it  is  rotating. 
Then  carrying  out  the  differentiation  in  (12)  using  (13)  and 
(11)  we  get : 

_  _ 

VXYZ  =  vxyz  +  a  I  (R  sin  $  +  y)  x  +  (R  cos  4>  -  x)y  (14) 

From  (14)  we  can  compute  Vxyz  (the  non-inertial  velocity) 
from  VXyz  (the  inertial  velocity)  and  the  dynamics  of  the 
airfoil.  In  the  above  derivation  we  could  also  have 
obtained : 

-  -  •  A  A 

RXYZ  =  =  Ro<sinj>x  +  cos<j>y) 

dtXyz 

Then : 

_  •  A  A  •  a 

RXyz  =  R[(ccsin4>  +a z  cos  <j)  )x  +  (a  cost})  -a  z  sint)))y  (16) 


Also: 


xy  z 


2  a  (v x  -  uy ) 

(17) 

a  (yx  -  xy) 

(18) 

yy) 

(19) 

,  we  get: 

Substituting  (16)  thru  (19)  in  (2),  we  get: 

fcf  T  dA  +  ///b  p  dvol  -  ///  {{R  (a  sin  <f>  +  a  ^  cos<}>)+2av+  a  y 

c.s.  c.v.  c.v. 

-  a  x}x  +  {  R (acos<}>  -  azsin<}))  -  2au  -  ax  -  a  y}  y  Jr p  dvol 


V  (p  V  •  dA) 
xyz  xyz 


+  a 


f  r  r 


(20) 
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Its  velocity  is  then,  by  definition: 

VXYZ  e  dP  (7) 

dtXYZ 

where  the  time  rate  of  change  is  that  seen  from  the  inertial 
reference  frame  (XYZ).  Then: 


VXYZ  =  dP  =  d  (R  +  r)  =  dR  +  dr 

dtXYZ  dtXYZ  dtXYZ  dtXYZ 

But  we  can  write,  from  figure  3: 


R 

=  -R  cos  <t>  x  +  R  s  i  n  <j>  y 

(9) 

r 

=  xx  +  yy 

(10) 

^xyz/XYZ 

A 

=  -  otz  (positive  c.c.w.) 

(ID 

Substituting  (9)  and  (10)  into  (8): 


VXYZ  =  d_  [-R  cos  <})  x  +  R  sin  <f>  y  +  xx  +  yy] 
dtXYZ 


=  d_  [  (-R  cos  <()  +  x)x  +  (R  sin<}>  +  y)  y] 

dtXYz  (12) 


As  the  quantity  in  bracket  is  not  expressed  in  the  same 
reference  frame  as  that  in  which  the  time  derivative  is 
taken,  we  need  to  use  the  general  rule  for  the  derivative  of 
a  vector: 


d  * 


dt 


XYZ 


xyz 


dt 


xyz 


-xyz/XYZ  x 


xyz 


xyz 


(13) 
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flf2 


xy 


reference  frame  is  a  non-inertial  reference 
frame  with  its  origin  at  the  rotation  point 
and  rotating  with  the  airfoil. 

reference  frame  is  a  non-inertial  reference 
frame  with  its  origin  fixed  at  a  point  on 
the  surface  of  the  airfoil  and  its  x-axis 
tangent  to  the  surface  while  the  y-axis  is 
perpendicular  to  the  surface. 


The  transformations  between  the  different  reference 


frames  are: 


and : 


/\  /s  ^ 

X  =  (cosa)  f i  +  (sinoc) 


V  =  (-sina)  f  ^  +  (cosa) 


A  A 

fj  =  (cosi|j)x  +  (-sin^)y 


(3) 


K  =  (sin^)x  +  (cos^)y 


(4) 


Substituting  (4)  into  (3)  and  using  the  trigonometric 
identities  cos(a-i|;)  =  cosa  cos  ij;  +  sina  sin  ij; 

sin  (a  -  ip)  =  sina  cos  -  cos  a  sin^ 


will  give: 


X  =  [cos  (a  -  ij; )  ]  x  +  [  s in  (  a  -  )  1  y 

Y=  [-sin(a-^)]x+  [cos ( a  -  ^  ] y 
Z  =  z 


(5] 


Consider  point  P  in  figure  3.  Its  position  vector  with 
respect  to  point  0  is: 


P  =  R  +  r 


(6) 
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Figure  3.  Dynamics  of  a  Pitching  Airfoil 
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In  these  two  basic  equations,  the  different  quantities 
have  the  following  meaning  (refer  to  figure  2): 


T  :  surface  forces 

dA  :  an  outward  normal  vector  to  the  surface 

element  of  magnitude  dA 

B  :  body  forces  in  an  inertial  c.v. 

Vxyz  :  velocity  vector  as  measured  from  the 

non-inertial  c.v. 

D/3  t  z  :  time  rate  of  change  as  measured  from  the 

w  non-inertial  c.v.  ~  ‘ 

dA  :  a  surface  element  of  the  control  surface 

R  :  position  vector  from  the  origin  of  the 

inertial  reference  frame  to  the  origin  of 
the  non-inertial  one 

oj  :  the  angular  velocity  of  the  non-inertial 

reference  frame  with  respect  to  the 
inertial  reference  frame 

r  :  position  vector  from  the  origin  of  the 

non-inertial  reference  frame  to  the  point 
under  consideration  in  the  non-inertial 
reference  frame 

We  will  now  apply  equations  (1)  and  (2)  to  a  specific 
control  volume  attached  to  an  airfoil  pitching  in  a  steady 
air  stream.  We  now  refer  to  figure  3. 

On  figure  3  we  note  that: 

XY  :  reference  frame  is  the  inertial  reference 
frame  with  its  origin  at  the  rotation 
point . 
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I I .  Theory  Development 


Introduction 

The  flow  about  a  pitching  airfoil  is  solved  analytically 
using  a  von  Karman-Pohlhausen  momentum  integral  technique. 
Because  our  control  volume  is  attached  to  the  airfoil,  the 
basic  equation  for  continuity  and  momentum  must  first  be 
derived  for  a  non-inertial  control  volume.  Finally,  the 
velocity  at  the  edge  of  the  boundary  layer  is  obtained  from 
the  inviscid  solution. 

Non-inertial  Control  Volume 

We  wish  to  develop  the  continuity  and  momentum  equations 
for  a  control  volume  attached  to  the  surface  of  an  airfoil 
undergoing  a  pitching  motion  in  a  steady  stream. 

For  the  most  general  control  volume,  the  continuity 
equation  is  given  by  (7:141-145): 


#  P  Vxvz  *  dA  =  -3 _  ///  P  dvol 

C.s.  *2  3txyz  c.v. 


(1) 


and  the  linear  momentum  equation  by: 


j tf  TdA  +///  Bp  dvol  -///[  R+2'u  x  V  +2to  x  rho  x  (w  x  r)]p  dvol 
c.s.  c.v.  c.v.  xyz 


■  tf  v(p  v*yz  dA)  +  —  /// Vv„ 7(pdvoi) 


C.S. 


3txyz  C-V< 


xyz 


(2) 
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he  investigated  the  pitching  airfoil  problem  and  empirically 
derived  rules  for  predicting  the  effect  of  a  constant  pitch 
rate  on  the  lift  of  airfoils  of  different  thicknesses  and 
camber  ratios. 

In  this  study,  the  viscous  flow  analysis  developed  by 
Lawrence  is  used;  however,  the  potential  flow  solution  makes 
use  of  a  new  concept  originating  from  a  jump  condition 
observed  by  Tupper.  The  rotational  motion  of  the  airfoil  is 
modeled  by  an  induced  camber.  Although  a  closer  match  to  the 
experimental  results  is  desireable,  the  thrust  of  this  study 
.  s  toward  achieving  better  physical  insight  into  the  problem 
rather  than  a  perfect  match  of  the  experimental  results.  To 
this  end,  some  assumptions  are  made  to  alleviate  the 
mathematical  difficulties  or  complexities  which  could  obscure 
the  problem. 


5 


Deekens.Kuebler,  Daley  Experiments  Data  Range 
Lawrence  Mass  Ingestion 

Lawrence  MRS 
.Pitching  Motion 


Results  of  Experimental  Investigation 
Stall  by  Kramer,  Deekens/ Kuebl er , 
Lawrence 


of  Dyna 
Daley 


(1).  In  his  experimental  set-up,  the  airfoil  was  fixed  in 
the  wind  tunnel  between  a  set  of  movable  guide  vanes;  he  used 
those  vanes  to  make  a  constant  change  of  the  angle  of  attack 
on  the  airfoil. 

In  1978,  Deekens  and  Kuebler  (2)  and  in  1982,  Daley  (3) 
conducted  experimental  studies  on  a  pitching  airfoil  in  a 
constant  velocity  smoke  tunnel.  Their  results  are  combined 
in  figure  1.  There  is  a  considerable  difference  between  the 
gust  experiment  and  the  pitching  airfoil  one.  Docken  (4)  has 
related  the  increase  in  the  stall  angle  of  attack  of  the  gust 
problem  to  a  reduction  of  the  adverse  pressure  gradient 
imposed  on  the  boundary-layer. 

In  1983,  Lawrence  (5)  investigated  in  detail  the 
pitching  airfoil  problem  using  a  von  Karman-Poh lhausen 
momentum-integral  technique  and  related  some  of  the  increase 
in  the  stall  angle  of  attack  to  the  same  reduction  of  the 
adverse  pressure  gradient;  further  portions  of  the  gain  were 
related  to  a  change  in  flow  separation  criterion  due  to  the 
moving  boundary  (Moore-Root-Sears  models)  and,  mostly,  mass 
introduction  into  the  boundary-layer  from  the  external  flow. 
His  results  are  also  included  in  figure  1. 

During  the  same  year,  Tupper  (6)  made  an  analytical 
study  of  the  potential  flow  about  an  airfoil  starting  from 
rest.  He  modeled  the  unsteadiness  of  the  flow  by  using 
vortex  shedding  at  the  trailing  edge.  Using  the  same  method, 


the  stall  occurs  at  a  significantly  greater  angle  of  attack 


and  is  much  more  severe  than  is  'ie  case  when  the  stall  angle 
of  attack  is  very  slowly  approached  (i.e.  static  stall 
conditions).  The  contributing  factors  to  this  phenomenon  are 
both  viscous  and  inviscid;  they  depend  on  pitch  rate,  free- 
stream  velocity,  airfoil  geometry,  Reynolds  number  and  Mach 
number  (10:304). 

Problem  Statement 

The  purpose  of  this  thesis  is  to  improve  the  knowledge 
of  the  physical  contributors  to  the  phenomenon  of  dynamic 
stall.  A  modified  von  Karman-Pohlhausen  momentum-integral 
method  is  used  to  analyze  the  flow  within  the  boundary  layer, 
while  a  complex  potential  flow  theory  is  used  to  analyze  the 
inviscid  flow.  Conformal  mapping  is  used  to  obtain  the  flow 
around  the  airfoil  from  that  of  a  circular  cylinder.  By 
using  these  relatively  simple  techniques,  it  is  easier  to 
keep  track  of  the  underlying  physics  of  the  problem. 

This  study  is  restricted  to  the  analysis  of  a  pitching 
airfoil  in  a  steady  flow  field;  as  the  control  volume  used  to 
study  the  boundary  layer  is  attached  to  this  pitching 
airfoil,  non-newtonian  fluid  mechanics  is  used. 

Background 

In  1932,  Max  Kramer  published  the  results  of  his 
experimental  study  of  a  dynamic  stall  induced  by  a  wind  gust 


INVESTIGATION  OF  POTENTIAL  AND  VISCOUS  FLOW 
EFFECTS  CONTRIBUTING  TO  DYNAMIC  STALL 

I .  Introduction 

Discussion 

The  phenomenon  of  dynamic  stall  has  received 
considerable  attention  from  fluid  dynamicists  since  1930. 
However,  most  of  the  research  has  been  experimental  and  the 
results  empirical.  Unsteady  aerodynamics  has  been  developed 
mainly  for  small  amplitude  oscillatory  displacements  about  an 
equilibrium  position.  Dynamic  stall  is  obviously  important 
to  the  helicopter  and  the  compressor  industries;  here  again, 
however,  the  amplitude  of  the  movements  about  an  equilibrium 
angle  of  attack  is  small  enough  to  allow  the  use  of 
perturbation  theory.  Consequently,  there  is  still  much  to  be 
learned  about  the  dynamic  stall  problem. 

In  this  study,  we  analyze  the  behavior  of  the  unsteady 
flow  pattern  about  an  airfoil  that  is  pitched  from  steady 
prestall  conditions  until  dynamic  stall  is  reached.  The 
total  angle  of  attack  excursion  may  be  as  large  as  40  to  50 
degrees,  making  any  use  of  perturbation  theories  impossible. 
When  the  airfoil  is  stalled  under  such  a  dynamic  condition. 


Dividing  thru  by  p<^x,  we  get: 

-  1  9£  dy  -  -  J*1  {  R(asincj>  +  a2cos<j>)  +  2av  +  ay  -  a2x  }  dy 

p  o  9x  p  o 

=  /h  9_  (uu)dy  -  ufi  fh  9u  dy  +  /h  9u  dy  (24) 

o  9x  o  9x  o  3t 

where  the  order  of  integration  and  differentiation  has  been 
interchanged.  We  note  that  in  (24)  all  time  derivatives  and 
velocities  are  those  measured  from  the  non-inertial  reference 
frame  (x,y)  . 

At  this  point,  we  need  to  derive  the  non-inertial  Euler 
equation  to  solve  for  9P/ 9x.  In  Appendix  A,  we  show  that 
this  equation,  when  applied  to  the  streamline  at  the  edge  of 
the  boundary  layer,  is  given  by: 


-  1  9P.  =  ue  9uP  +  9ug  +  {  R(dsint))  +  a2cos<f>)  +  2a ve  +  aS  -  a2x}  (25) 
p  9x  9x  9t 

Substituting  (25)  into  (24)  and  grouping  terms,  we  get: 


^-/h  [ue  iMs.  +  9_(ue  "  u)  +  2a(v_  -  v)  +  a(6  -  y)]  dy 
p  o  9x  9t 

=  -/h  9_  (uu)dy  +  ue  /h  9u  dy 
o  9x  o  9x 


It  can  easily  be  shown  that  (5:64-65) 


and 


/h  3_  (ue  -  u)dy  =  a__  (u e6i) 

o  3t  at 

where  6-^  and  62  are  the  displacement  and  momentum  thickness, 
respectively. 

Therefore  (26)  can  be  written  as: 

TW  =  i_  (ue2<S2)  +  ue  /  3Ua  U  +  iL  (“e  $1  )  +/h[2a(ve-v)+a(6-y)]dy,27) 

p  ax  vax  /0l  at  o 

In  the  derivation  of  (27)  ,  the  only  assumptions  made  were: 

a)  rigid  airfoil; 

b)  incompressible  flow; 

c)  center  of  rotation  located  on  airfoil  chord; 

d)  two  dimensional  flow; 

e)  ap/ax  =  dP/dx;  and; 

f)  Weight  of  fluid  in  c.v.  negligible. 

Before  proceeding  to  transform  (27)  into  the  von  Karman- 
Pohlhausen  working  equation,  we  analyze  its  last  term: 

I  =  Jh  {2a  (ve  -v  )  +  a(6  -  y)}  dy  (28) 

0 

The  second  term  of  I  gives: 

f^a  (6-y)=a(6h-h2) 

0  2~ 

There  remains  tc  integrate  the  first  term.  As  a  first 
approximation,  we  can  treat  the  y-component  of  the  velocity 
as  a  Hiemenz  flow  (stagnation  flow  in  plane)  for  which  an 
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exact  solution  exists  (8:95-99).  The  velocity  distribution 


flow  in  the  vicinity  of 
stagnation  point  is 
given  by: 

u  =  a  x  Voo 
v  =-a  y  vro 


Figure  6.  Hiemenz  Flow 


Now 

2a  ( ve  -  v  )  dy  =  2a  ve  /^ ( 1  -_v_)  dy 
0  0  ve 


From  ( 30 ) : 


v  =  -a  y 

vro 


and 


6 


Therefore  v  _ 

ve  6 

So  that 


/h  2a (v  -  v )  dy  =2a  v  (h  -h^) 
o  25 


Therefore 


I  =  a(6h  -  h2)  +  2a(h  -  h2)ve 


the 


(30) 


(31) 


As  I  will  eventually  be  evaluated  at  y 
immediately : 


I  «  d  fi2  •  t 
|  +  a5ve 


=  6 ,  we  pose 


so  that  (27)  becomes: 

TW  3(u|  62)  +  uP(9up)6i  +  3_  (ue<Si)  +  a  62+  a6ve  (32) 

p  3x  3x  3t  2 

Development  of  the  Von  Karman-Pohlhausen  Method 

In  the  development  of  the  von  Karman-Pohlhausen  working 

equation,  we  follow  the  procedure  described  by  Schlichting 

(8 : 206-223) . 

If  we  expand  the  derivatives  in  the  momentum  integral 
equation  (27),  multiply  it  by  62/  vUe  and  use  the  cl°sure 
equation  (5:13-25),  we  get: 


62  ue  362  +  (2  +  \  6f  3Ue  +  1  61  6§  1  3ue 

v  3x~  \  6T /  v  3x  2  67  v  ue  9t 


+  I  52 
v  ue 

Adding  and  subtracting  the  quantity 


(2  +  1  3j  )  62  / 1_  3Ug\  +  (2  +  \  __6|_  {R (asi  n4>+a2cos4>) 

'  2  Ci 2  v  \  Up  3 1  /  \  62  /  v 


+  2ave  +  a6  } 


gives 


12  +  6j_\ 

6|  Eul  +  ue  62  362  + 

62 1  - 

(2  +  1 

6i\ 

§1 

M 

l  62  > 

v  v  3x 

vue 

\  2 

62/ 

V 

\ue  at/ 

=  LW  62 
M  Ue 
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where : 


Eul  =  R  (cxsi  n4>  +  a2cos<}>)  +  2ave  +  a6  +  Blip  +  1  Bug 


Bx  ue  3t 


We  now  assume  that  the  velocity  profile  within  the 
boundary-layer  has  the  form: 


u  =  A  +  Bn  +  Cn2  +  Dn3  +  En4 
ue 

where  g  =  y/6 /  and  is  subjected  to  the  boundary  conditions: 


at  y=0  :  u  =  0 
at  y=6  :  u  =  ue 


v32u  _  -(Eul )up 

~§yr" 

_3u  _  0  32l 

3y  By2 


We  define 


6 2  Eul 


and  solve  (37)  for  its  coefficient  by  applying  the  five 
boundary  conditions;  we  get 


The  velocity  profile  is  then  given  in  terms  of  A  by: 


u  =  (2n  -  2n3  +n4)  +  A  (n  -  3n2  +  3n3  -n4)  (41) 

ue  6 

By  applying  the  definition  of  the  displacement  and  momentum 
thickness,  we  get  (8:209) 


6 !  _  3  -  A 

6  10  120 

62  _  37  -  A  -  A2 

6  “  315  945  9072 


The  shear  stress  at  the  wall  is  given  by 

TW  =p3u 

9y  w=n 


so  that 


2  +  _A_ 

6 


If  we  assume  that  separation  occurs  when  the  shear  stress  at 
the  wall  is  zero,  then  from  (45),  A=  -12  at  separation. 
However,  Moore,  Root  and  Sears  have  developed  a  better  model 
for  separation  for  the  case  of  a  moving  wall  (13:123-126). 
In  Appendix  B,  we  show  that  the  value  of  j\  at  separation  is 
given  by  the  solution  of: 


(As  +12)2  ( 9A 2  -  136A  +  528) 

5 1 2 ( As - 6 ) 3 


“WALL 


where  u 


WALL 


is  the  tangential  component  of  the  wall  velocity 


at  the  separation  point. 
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Let  us  now  define: 


Z  Eul 


Then,  from  (39),  (43)  and  (48)  the  shape  factors  K  and  A 

satisfy  the  universal  relation: 


37  -  A  -  A  2  \  A 

315  945  9072  ) 


Let  us  also  define 


fl (K)  E  6j_  =  /_3  -  _A_\  /_37  -  JV 


10  120 )  \  315  945  9072 


y  ue 


2  +  A  \  /_37  -  _A_  -  A  2 
6  /  V  31 5  945  9072 


F  (K)  =  2f2(K)  -  4K  -  2Kfx(K) 


and  note  that  (from  (47)) 


Up  £2.  Mi  =  1  Up  dZ 
v  3x  2  dx 


Then,  if  we  substitute  (47)  thru  (53)  in(35)  and 

rearrange  some  terms,  we  get: 

dZ  =  _1^_  {  F(K)  +  {4  +  fi(K)}  _Z  btia  +  R (cisi n4>  +  a2coS(j>) 


dx  ue 


ue  t 

+  2ave  +  a6  -  262!  } 

ue 


Z  =  K(Eul)"1 


where  again,  Eul  is  given  by  (36).  In  this  set  of  equations, 

the  quantities  a,  a,  R  and  cj>  are  known  from  the  geometry  and 

the  dynamics  of  the  problem;  the  quantities  F(K),  f^K), 

f2  (K)  are  functions  only  of  the  shape  parameter  K  (or  A  )  ; 
finally  the  other  quantities  ue,  ve,  dU e/dx,  dUe/dt  ,  I  and 
Eul  are  derivable  from  the  inviscid  flow  outside  the  boundary 
layer.  However,  two  problems  remain  to  be  solved.  First,  to 
evaluate  I,  we  need  to  make  some  assumption,  such  as  the 
Hiementz  flow  assumption  of  page  17.  In  this  case: 

26 2 1  =2  6l  6  / ac+  ave\  =  2  Zf3(K)  /a6  +  avd  (56) 

u  Ue  ue  u  62  \ 2  /  ue  V  2  / 

where  f3(K)  =  6_  (57) 

62 

Second,  5  occurs  directly  in  the  term  a6  of  (54)  and  in  Eul. 
We  therefore  have  a  closure  problem  as  the  value  of  6  is 
required  not  only  for  these  terms  but  also  to  obtain  the 
value  of  Ve  from  inviscid  flow.  A  computer  iteration  routine 
could  be  used  to  evaluate  <5(x).  As  we  start  from  steady 
state  conditions  before  subjecting  the  airfoil  to  a  pitching 
motion,  a  = a  =0  and  we  recover  the  steady,  non-inertial 
von  Karman-Pohlhausen  working  equations  from  (54)  and  (55). 
Their  integration  will  provide  the  initial  boundary- layer 
thickness  distribution  6(x).  As  the  pitching  motion  start, 
(54)  and  (55)  are  solved  using  the  previous  S(x)  and  so  on 
till  separation  occurs  at  the  quarter  chord.  This  should 
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yield  an  exact  solution  of  the  pitching  airfoil  dynamic  stall 
problem  within  the  assumptions  made  on  pages  16  and  17  .  An 
offset  of  this  approach  would  be  the  study  of  the  change  of 
the  boundary-layer  thickness  due  to  the  pitching  motion,  and 
the  effect  of  this  change  on  the  dynamic  stall. 

An  alternative  approach  is  that  of  Lawrence  (5).  He 
showed  that  the  effect  of  the  hypothetical  (non-iner t ia 1 ) 
body  forces  on  a  control  volume  where  the  upper  limit  is  the 
edge  of  the  boundary  layer  :an  be  neglected.  Then  he  used 
the  unsteady  but  inertial  Euler  equation: 

ue*au£  +  =  i  dP  (58) 

3x  3t  p  dx 

where  Ue* is  the  velocity  viewed  from  the  inertial  reference 
frame,  and  assumed  that 

ue  =  ue*  -  Rasin<{>  (59) 

He  finally  accounted  for  the  non-inertial  conditions  by 
assuming  that  mass  was  ingested  in  the  boundary  layer  through 
the  superposition  of  a  Hiementz  flow.  His  development  of  the 
Von  Karman-Pohlhausen  method  leads  to  the  equations: 

dZ  ={  F(K)  +  [4  +  f i ( K )] _ Z  cHJg_  +  [4  +  2fx(K)  -  2f 3 (K)]  • 

dx  ue  at 

ZRasimj)  3uP  +  2/c  \zf3(K)  Racos<{>}  (60) 

ue  9x  \  5  /  ue 

and 

Z  =  K  ["(Rasing)  1  aue  +  3uQ  +  1  auel 


(61) 


V w.  VT  v„  V  V  V  V  ' 


For  the  Joukowski  pitching  airfoil  and  the  Hiementz  flow 
model,  he  used  C  =  4  6  .  His  model  achieved  considerable 
improvement  in  the  insight  of  factors  contributing  to  the 
dynamic  stall  effects.  As  it  has  the  advantage  of  not 
requiring  any  a  priori  knowledge  of  6(x)  ,  it  is  used  in  the 
rest  of  this  study. 

Inviscid  Solution 

In  order  to  solve  the  modified  von  Karman-Poh lhausen 
equations,  we  need  to  evaluate  the  velocity  at  the  edge  of 
the  boundary- 1 ayer ,  Ue.  In  this  section,  we  used  potential 
flow  theory  and  conformal  mapping  to  evaluate  it. 

Figure  7  represents  the  flow  of  a  circular  cylinder  of 
unit  radius  undergoing  rigid  body  rotation.  The  vorticity  is 

given  by: 

w  =  curl  V  (62) 

For  rigid  body  rotation: 

Ur  =  0  U0  =  -  a  r  (63) 

Figure  7.  Rigid  Body  Rotation 

The  minus  sign  occurs  because  0  is  positive  counter¬ 
clockwise,  therefore 

a  =  -2a  k  (64) 
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For  any  given  time,  v  and  v'  are  constant  and  tangent  to 
their  respective  circles.  Applying  Stokes's  circulation 
theorem  to  these  two  circles  gives,  and  as  ds  =  r  d  0  : 

2 

9  v'  ds  =  -2am  r‘  r'  <  a  (65a) 

r' 


and 


0  v  ds  =  -2a,7T  r  >  a  (66a) 

r 


The  second  integral  is  evaluated  at  r  =  a  =  1  as  only  the 
cylinder  is  rotating  and  therefore  the  area  between  r  =  1  and 
r  =  r  contributes  nothing  to  the  integral.  Since  v'  and  v 
are  constant  on  their  respective  circle,  we  also  have 

9  v1  ds  =  v‘  2-nr' 
r' 

9  v  ds  =  v  2nr 
r 

By  comparing  (65)  and  (66)  we  find: 

v  1  =  -  a  r  ' 

we  have  a  rigid  rotator  for  r  <  a;  and 

v  =  -  a  /r 

a  vortex  centered  at  r  =  0  and  of  strength  rR  -  +2  ma  when 
r  >  a. 

We  conclude  that  the  flow  outside  a  rotating  cylinder 
may  be  accurately  modelled  by  superposition  of  the  flow  about 


(65b) 

(66b) 
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a  non-rotating  cylinder  and  the  flow  induced  by  a  vortex  of 
strength  +2ira  located  at  the  center  of  the  circle.  In  order 
to  satisfy  Helmoltz's  vortex  Laws  (17:54),  we  must  assume 
that  a  vortex  of  equal  and  opposite  strength  simultaneously 
appears  at  infinity. 

The  stream  function  (  ip  )  and  the  potential  function  (  <j>  ) 
for  the  flow  about  a  non-rotating  cylinder  with  circulation 
in  a  free  stream  of  uniform  velocity  and  having  an  angle 
of  attack  a  are  in  polar  coordinates  (r,0)  given  by  (12:84): 


=  r  Uoo  sin(9  -  a)  /I  -  a2\  +  r  In  /  r 


\  +  r_  In  l  r 

I  2tt  la. 


4>  =  r  Uoo  cos(0  -  a)  /I  +  a2  \  -  T  0 

l  p"2"  / 

The  velocities  at  any  point  (r,6)  are  then 


Ur  =  Uoo  cos(0  -a  )  fl  -  a: 


ie  =  -  Uro  sin(9  -  a)  ( 1  +  a2 )-  r 

'  r 2 '  2TTr 


On  the  surface  of  the  cylinder,  r  =  a  and  therefore  Ur  = 
everywhere  (from  (69))  and: 


U9  =  -2Uoo  sin(0  -&  )  -  r 

2Tia 


We  apply  Lhe  Kutta  condition  in  order  to  have  a  stagnation 
point  at  the  trailing  edge;  i.e.  for  0=0,  we  obtain  the 
steady  state  vortex  strength: 


r 


=  4:'aU,.,  sinuss 


(71) 


1  ss 

We  now  superpose  the  vortex  due  to  the  rotation  to  the 
above  non-rotating  cylinder  in  steady  state  conditions.  The 
stream  and  potential  functions  become,  for  a  cylinder  of  unit 
radius: 


=  rlloc. 

s i n ( 6  -  a) 

I1'*) 

+ 

|  F  +  a  jin  r 

(72) 

4)  =  rU„, 

cos(0  -  a) 

/I  +  1  ' 

l  , 

i  ■  ( 

'  r  +  a 

>  2 "  > 

(73) 

Froiri  (72)  we  find  that  Ur  is  still  given  by  (69)  and 
therefore  remains  zero  everywhere  on  the  surface.  But  now: 


Ue  =  -U„  sin(0  -  ex)  /I  +  1  \  -  /_£_  +  i\  l 

\  r2 1  l  2tt  /  r 


In  this  study,  group  together  the  two  circulation  terms 

so  that  JT  =  _r_  +  a  ,  and  TgS  is  still  given  by  (71).  The 

2n  2  n 

reason  we  choose  to  group  these  terms,  is  to  obtain  results 
which  can  be  compared  with  those  of  Tupper. 


As  rotation  occurs,  the  angle  of  attack  is  continuously 
changing.  From  (71)  we  see  that  this  continuously  increases 
Fss  ;  this  particular  augmentation  does  not  occur 
instantaneously,  however.  Tupper  (6;5-18)  has  developed  a 
procedure  of  vortex  shedding  at  the  trailing  edge  to  account 
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for  the  progressive  growth  of  the  circulation  following  an 
impulsive  start  of  angular  rotation. 

In  this  method,  a  vortex  is  shed  in  the  wake  of  the 
airfoil  behind  the  trailing  edge.  At  the  same  time  a  bound 
vortex  of  equal  and  opposite  strength  is  inserted  inside  the 
airfoil.  The  location  of  the  shedded  vortex  is  U»  At  where 
At  is  the  time  elapsed  since  the  start  of  the  impulsive 
motion;  the  bound  vortex  is  located  inside  the  airfoil  in 
such  a  way  that  the  surface  of  the  airfoil  remains  a 
streamline  and  that  the  Kutta  condition  is  maintained. 


In  terms  of  the  flow  about  the  circular  cylinder,  which 
will  eventually  be  mapped  into  the  airfoil,  it  can  be  shown 
that  the  Kutta  condition  needs  to  be  applied  at  0  =  0  (11  :469) 
and  that  the  location  of  the  pair  of  vortices  is  obtained 
from  the  circle  theorem  (12:84-85).  This  way,  referring  to 

Figure  8,  if  the  shedded  vortex 
y'  is  l°cated  at  (r^,0^),  the 

/  surface  remains  a  streamline  if 


Figure  8.  Circle  Theorem 


the  bounded  member  of  the  pair 
is  at  ( 1  / rj  ,  6j  )  ■  The 
contribution  of  this  pair  to 
the  stream  function  is  then: 


In  [_1  1  +  r.j2r*  -  2n 

L r i 2  rl  +  r-j z  -  2r-j 


r  cos(9-  9- 
rcos(0  -  0i , 
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phenomenon,  then  the  relaxation  coefficient  used  in  the  mass 
ingestion  model  is  not  four  but  something  larger.  Perhaps 
this  relaxation  coefficient  should  be  a  function  of  position 
along  the  airfoil.  Larger  coefficients  could  be  rationalized 
near  the  leading  edge  where  relaxation  distances  are  smaller. 

In  this  study,  the  induced  camber  was  assumed  to  build 
up  instantaneously  to  its  value  specified  by  the  non- 
dimensional  pitch  rate;  however,  pitch  rate  does  take  a 
finite  time  period  to  attain  its  steady-state  value.  This 
effect  would  certainly  reduce  the  reported  induced  camber 
effect,  although  the  magnitude  of  the  reduction  could  only  be 
estimated  from  more  precise  knowledge  of  the  experimental 
apparatus  used  by  Deekens  and  Kuebler  (2). 

Thickness  Effect 

The  mass  ingestion  model  program  was  run  for  symmetrical 
Joukowski  airfoils  of  different  thickness  ratios,  represented 
by  the  value  the  real  part  of  in  the  plane  transformation 
given  by  equation  (82).  The  effect  of  thickness  on  the 
increase  in  separation  angle  of  attack  (A<*sep)  due  to  the 
pitching  motion  is  shown  in  figure  1  3.  Ac*sep  increases  with 
thickness  for  all  pitch  rates. 


For  thin  airfoils  (real  part  of  m  on  the  order  of  0.05 
or  less),  separation  occurred  at  the  leading  edge  under 
static  conditions.  The  dynamic  effect  could  therefore  not  be 
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— -  Mass  Ingestion 

- - Induced  Camber 

a:  per  Lawrence 
b:  with  correction  for  Kutta 
condit  ion 
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Figure  12 


Induced  Camber  Model  Resul 


Our  induced  camber  model,  is  then  applied  as  follows: 

a)  as  the  model  represents  the  rotating  motion  of 
the  airfoil,  it  replaces  the  mass  ingestion  model 
used  by  Lawrence  (5:47-51);  its  effect  is  therefore 
added  to  the  MRS  line  of  figure  1; 

b)  from  Appendix  C,  the  camber  induced  by  the 

rotating  motion  of  an  airfoil  of  type  J015  is: 

f/c  -  3  -4  7  ot 

therefore,  for  any  given  oe^,  the  increase  in 
static  stall  angle  of  attack  is  obtained  from  the 
data  of  figure  10  and  added  to  the  MRS  curve  for 
the  same  oND> 

The  result  is  shown  in  figure  12.  The  induced  camber 
effect  is  very  large,  much  larger  than  expected.  A  review  of 
the  data  compiled  for  laminar  airfoils  by  Althaus  (18:13-20) 
does  not  show  such  an  increase;  the  fact  that  a  Joukowski 
airfoil  has  a  circular  arc  profile  with  its  maximum  close  to 
the  half  chord,  while  laminar  profiles  has  parabolic  camber 
with  maximum  close  to  the  quarter-chord  could  explain  this 
discrepency.  It  was  anticipated  that  the  induced  camber 
effect  would  fall  short  of  the  experimental  data  range  when 
superimposed  on  the  pitch  and  MRS  effects.  If,  as  we  think, 
induced  camber  and  mass  ingestion  model  the  same  physical 
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Camber  Effect 


As  the  results  of  Lawrence  (5:53)  fall  short  of  the 
experimental  data,  an  effort  was  made  to  find  other 
contributing  factors  to  the  dynamic  stall.  Tupper's  results 
(6:56)  for  a  rotating  cylinder  showed  an  unexplained  jump  in 
lift  coefficient  immediately  after  the  start  of  pitching 
motion;  as  a  camber  airfoil  has  such  a  jump  over  its 
corresponding  symetrical  airfoil,  it  was  hypothesized  that 
the  rotating  motion  could  be  modeled  by  an  induced  camber,  as 
shown  in  Appendix  C. 

The  program  was  then  run  for  airfoil  of  different  camber 
ratio  but  having  the  same  thickness  ratio.  The  result  was 
that  camber  had  no  noticeable  effect  on  the  curves  Aasep 
versus  ca/2U„for  both  the  MRS  and  mass  injection  cases  shown  in 
figure  1.  However,  the  static  quarter-chord  separation  angle 
of  attack  sharply  increases  with  camber  ratio  as  shown  in 
figure  10,  where  the  camber  ratios  were  measured  from  the 
plots  of  the  airfoil  shown  in  figure  11  as  a  function  of  the 
vertical  displacement  of  the  mapping  circle,  i.e.  the 
imaginary  part  of  the  complex  number  m  of  the  plane 
transformation  given  by  equation  (82). 

To  obtain  the  true  static  separation  angle  of  attack  for 
a  cambered  airfoil,  the  value  of  asep  obtained  from  the 
program  had  to  be  reduced  by  the  value  of  p  used  in  the  plane 
transformation  given  by  equation  (82);  this  difference  is 
strictly  due  to  the  geometry  of  the  transformation. 
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III. 


Results 


Preliminaries 

Before  attempting  to  progress  from  the  results  of 
Lawrence  and  Tupper,  several  modifications  were  made  on 
Lawrence  program  (5:91-96).  Firstly,  it  was  generalized  to 
compute  the  boundary  layer  of  a  general  Joukowski  airfoil  of 
arbitrary  thickness  and  camber  ratios;  secondly,  it  was 
modified  to  accept  any  location  for  the  point  of  rotation  of 
the  airfoil;  thirdly,  a  procedure  was  incorporated  to  compute 
the  quarter-chord  separation  criterion  in  accordance  with  the 
MRS  model  developed  in  Appendix  B;  and  finally,  a  correction 
had  to  be  made  in  the  subroutine  computing  the  velocity  to 
properly  account  for  the  Kutta  condition  at  the  trailing 
edge. 

The  program  was  then  run  for  a  symmetric  Joukowski 
airfoil  of  thickness  ratio  0.15  at  zero  degree  angle  of 
attack.  The  results  matched  those  published  by  Schlichting 
(8:213)  for  the  same  test  conditions,  providing  assurance 
that  the  program  was  working  properly. 

Different  combinations  of  thickness  and  camber  could 
then  be  studied  as  a  function  of  a  non-dimensional  pitch 
parameter  (0.5  cd/2U«) .  It  has  been  shown  that  when  this 
parameter  is  used,  the  individual  values  given  to  c,  6  or 
Uw  do  not  effect  the  results  as  long  as  they  combine  into  the 
same  parameter  value. 
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The  required  values  ofdUe/3s  and  Ue  where  obtained  from  the 
non-inertial  value  of  the  first  part  and  using 


where  A0  =  2  and  As  =  ( (y (0  +.01))  -  y(e-.Ol))2  +  (x(0  +.01)  -  x(e  -0l)2)^ 
For  this  part,  the  integration  was  started  with  the  airfoil 
at  any  initial  angle  cr  attack  without  pitching  motion  to 
allow  the  flow  to  stabilize  around  the  forward  stagnation 
point.  The  pitching  motion  was  then  started,  and  the  von 
Ka r ma n- Poh 1  ha us en  equations  were  integrated  from  the 
stagnation  point  to  the  quarter-chord.  The  MRS  separation 
criterion  developed  in  Appendix  B  was  used  to  check  for 
separation  at  the  quarter-chord;  if  separation  did  not  occur, 
the  solution  was  repeated  using  a  larger  initial  angle  of 


attack . 


of  attack;  this  way,  for  a  symmetrical  airfoil,  the  steady- 
state  circulation  is  zero,  and  steady-state  conditions  exist 
from  the  start.  The  velocity  at  the  surface  of  the  airfoil 
was  computed  from  equations  (70)  and  (71)  for  any  required 
angular  position  in  the  plane.  3u/3e  was  also  computed  using 
finite  differences  as  3u/36  =  (U  (  e  +  .01)  -U  ( 8 -.01) ) /Ae  , 
where  Ae =  .02.  For  this  first  step,  9u/3t  is  zero 

everywhere  on  the  surface,  i.e.  steady-state  condition.  The 
airfoil  was  then  continuously  pitched  by  increments  Aa = 
a  At  =  0.01  degrees.  For  each  a  ,  the  velocity  U,  and  the 
velocity  gradient  3u/3e  were  computed  as  for  the  case  of 
a =  0.  The  local  time  rate  of  change  of  velocity  was  also 
computed  by  finite  differences  as 

9u  _  u(«  +  Aa )  -  u(a  -  Aa) 

3t  '  2At 

where  At  =  Aa/a.  A  subroutine  was  then  used  to  compute  the 
non-inertial  velocity  (using  equation  (14))  from  the  above 
inertial  velocity. 

For  the  second  part,  the  von  Karman-Pohlhausen  working 
equations  (59)  and  (60)  where  numerically  integrated  using 
the  isoclines  method.  The  initial  conditions  required  are 
(8:211) : 

K0  =  0.0770 
ZD  =  0.0770/ (  u/  x ) 0 
(dZ/dx) Q  =  -0.0652  (  2U/  x2)0/(  u/  x)2 


I 
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The  velocity  in  the  Z  plane  will  be 


dF  =  dF  d£  d£* 
dz  dp  dp'dz 


but  from  (82) 


de.'=  e-iB 

dp  e 

which  has  a  magnitude  equal  to  1;  and  from  (83) 


dz  _  1  -  /  pt1  \ 
dp'  ~  [  p*  j 


The  velocity  distribution  in  the  Z  plane  is  obtained  by 
substituting  (87)  and  (88)  into  (86),  to  get: 

w ( z )  =  e+lB  1  -  1  w(p )  (89) 

where  w(  p)  =  u(p  )  -  i  v(  p)  is  known  from  equation  (77)  and 
(78).  As  we  are  only  interested  in  the  velocity  on  the 
surface  of  the  airfoil  in  the  Z  plane,  only  the  magnitude  of 
w(z)  is  required;  its  direction  must  be  tangent  to  the 
surface  everywhere. 

We  now  have  all  the  elements  required  to  solve  the 
pitching  airfoil  problem. 


Numerical  Solution  Process 


The  numerical  solution  to  the  pitching  airfoil  problem 
was  done  in  two  parts.  First  the  velocity  distribution  at 
the  surface  of  the  airfoil  was  obtained.  We  started  with  the 
non  rotating  airfoil  sitting  in  a  steady  flow  at  zero  angle 
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p  Plane 


p'  Plane 


Z  Plane 


center  is  at  the  origin  of  the  p  plane.  If  a  point  in  this 
plane  is  given  by  p  =  r  e1  ,  then  it  can  be  transformed  to  a 
p'  plane  by  (11:461) 


The  p'  plane  has  all  the  points  of  the  p  plane  rotated  g 
radians  clockwise  and  displaced  a  distance  y.  It  therefore 
represents  the  flow  about  a  rotating  cylinder  with  its  center 
located  at  y,  a  complex  number.  The  flow  may  now  be 
transformed  to  that  about  a  Joukowski  airfoil  (the  Z  plane) 
by  the  transformation. 


Z  =  p'  +  Pt 


The  point  is  that  point  of  the  p'  plane  where  the 

cylinder  circumference  crosses  the  positive  x-axis;  it  maps 
to  the  trailing  edge  of  the  Joukowski  airfoil. 

Knowing  the  velocity  distribution  in  the  p  plane  and  the 
mapping  functions  from  the  p  to  the  pr  and  on  to  the  Z  plane, 
we  can  find  the  velocity  distribution  in  the  Z  plane.  If  <J> 
and  ijj  are  the  potential  and  stream  functions  associated  with 
the  flow  in  the  p  plane,  we  can  form  the  complex  potential 


F(p )  =  <P  +  iip 


where  P  =  x  +  iy,  in  the  p  plane;  the  complex  velocity  in  the 


p  plane  is  then: 


w(p)  =  dF 


For  the  first  pair,  N  =  1  and  on  the  surface  r  =  1;  then 
applying  the  Kutta  condition  for  0=0  and  solving  for  >  we 
get : 


P i  =  4-itUoo  sin(a  -  since..)  (ri  -  1) 

ri  (rx  +  1] 


this  shedded  vortex  is  located  at 


r2  =  1  +  ILAt 

0i  =  0 


The  Kth  vortex  pairs  are  evaluated  similarly  once  the  first 
K-l  pairs  are  known.  In  general. 


rk  =  1  /rk  -  1  \  4irUoo  (sina  -  sinass) 

rk  \rk  +  *  /  k-l  x  (8i  \ 

-l  riV/rf-J,  ) 

i=l  \l+r-p  -  2rico s0 -j / 

Again  (rk,  q  ^ )  is  given  by  (80)  while  the  position  of  the 
previously  shed  vortices  is  given  by: 

ri  ’  *"i ,  t-At  +  urAt 

0i  =  9i,t-at  +  u0At 

where  Ur  and  are  computed  from  (77)  and  (78)  and 

indicates  position  before  the  time  increment. 


Conformal  Mappinq  -  Joukowski  Transf ormation 


Equations  (77)  and  (78)  give  us  the  velocity  at  any 
point  of  the  flow  field  about  a  rotating  cylinder  whose 
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After  a  second  time  increment  A  t,  a  second  vortex  is  shed 
at  Uoo  At,  while  the  first  vortex  has  moved  in  the  flow  by  a 
distance  uAt,  where  U  is  the  velocity  of  the  flow  at  the 
original  vortex  location  (r^,e^).  And  so  on  with  each  time 
increment.  At  each  step,  r can  be  computed  by  applying  the 
Kutta  condition  at  0=  0. 

The  complete  model  for  the  stream  function  of  our 
rotating  cylinder,  where  the  rotation  starts  ins tanteous  ly 
from  steady-state  conditions  of  Uoo  and  a  =  ass,  is  then: 

/  x  N 

ip  =  rUoo  sin(0  -  a)  (1  -  1  )  +  2Uco  sina  In  r  +[  i  (76) 

\  “r2/  i  =  l  ' 

where  the  ipr  are  given  by  (75).  As  time  passes,  the  strength 

'i 

of  the  vortices  remains  constant,  but  their  location  changes 
in  accordance  with  the  flow  velocity  at  their  prior  location 
in  the  field.  The  strength  of  the  first  pair  is  evaluated 
from  the  Kutta  condition.  Then  the  strength  of  the  following 
pair  may  be  computed  as  that  of  all  previous  pairs  is  known. 
From  (75)  and  (76),  we  get: 

Ur  =  9^  _  Uoo  cos(0  -  a)  /  1  -  1  \  (77) 

r  30'  l  ~rr  I 

so  that  for  r  =  1,  Ur  =  0  and  the  surface  is  a  streamline  as 

required.  Also: 

u0  =  -9^  _  -um  sin(0  -  a)  /I  +  1  \  -  ZU„sina<;<; 

9r  '  r2  /  r2 


N 

"  l  lili2  fri2  r  -  r-j  cos (9  -  8j )  -  r  -  rj  cos(0  -  9  - ) 

1=1  2tt  1  +  r2r^  -  2rricos(0  -  0i )  r2+ri  -2rricos (0-0-j ) 


(78) 
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studied.  This  result  is  compatible  with  the  well  known 
leading  edge  bubble  separation  (9:220). 

Location  of  the  Point  of  Rotation 

Another  factor  which  could  affect  the  results,  is  the 
location  about  which  the  airfoil  undergoes  its  pitching 
motion.  As  the  point  of  rotation  is  moved  toward  the 
trailing  edge,  the  stream  velocity  induced  by  the  pitching 
motion  is  increased,  for  the  same  value  of  c*Nq  when  compared 
to  rotation  about  the  mid-chord  point.  It  was  therefore 
expected  that,  on  one  side,  more  mass  would  be  ingested  at 
corresponding  points  on  upper  surface  of  the  airfoil  ahead  of 
the  rotation  point  and,  on  the  other  side,  the  MRS 
contribution  would  be  increased  as  the  wall  velocity 
increases.  Both  effects  would  tend  to  increase  the  quarter- 
chord  separation  angle. 

The  program  was  run  for  several  locations  of  the 
rotation  point  and  the  results,  shown  in  figure  14,  confirmed 
our  expectations. 
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PARAMETER  :  Point  of  Rotation  (jn  %  chord) 

A:  Portion  of  total  effect  due  to  mass  ingestion  (foraND  0.04) 


l.oc 

0.75c 

0.50c 

0.25c 


C  d  /  2  U* 


Figure  14.  Increase  of  Quarter-Chord  Separation 
Angle  of  Attack  With  Location  of  Rotation  Point 
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IV.  Conclusion  and  Recommendations 


Conclusion 

It  has  been  shown  that  as  an  airfoil  pitches  at  a 
constant  a,  the  separation  angle  of  attack  is  significantly 
increased.  This  increment  is  more  pronounced  for  thicker 
airfoils;  as  the  point  of  rotation  moves  from  the  leading  to 
the  trailing  edge  of  the  airfoil,  Actsep  also  increases. 

Camber  was  shown  to  have  no  significant  effect  on  the 
viscous  solution,  i.e.  the  value  of  AasepWhen  computed  as  the 
difference  between  the  dynamic  and  static  separation  angles 
of  attack  for  the  same  cambered  airfoil.  However,  a  very 
large  increase  in  static  stall  angle  of  attack,  when  compared 
to  that  of  the  same  thickness  but  uncambered  airfoil,  was 
obtained. 

An  effort  was  made  to  model  the  pitching  motion  of  a 
symmetrical  airfoil  by  an  induced  camber  due  to  rotation. 
The  results  obtained  show  a  substantial  increase  in  Aasep. 
While  the  mass  ingestion  model  results  are  below  the 
experimental  results,  the  induced  camber  model  gives  larger 
results.  If  the  hypothesis  of  induced  camber  is  correct, 
then  some  other  factors  contributing  to  the  dynamic  quarter- 
chord  separation  must  therefore  exist,  to  reduce  the  camber 
effect  to  the  experimental  data  range.  Camber  build-up  time 
delay  could  be  such  a  factor.  Separation  aft  of  the  quarter- 
chord  would  also  change  the  pressure  distribution  over  the 


airfoil;  in  this  study,  such  aft  separation  effects  have  been 
neglected  but  this  too  might  partially  explain  differences. 

Recommendations 

This  study  has  raised  many  unanswered  questions;  it  is 
recommended  that  the  following  be  further  studied: 

a)  A  computer  program  which  would  use  the  velocity 
distribution  about  a  circular  cylinder  shedding 
discrete  vortices  at  its  point  mapping  to  the 
trailing  edge,  as  developed  by  Tupper,  would  be 
more  representative  of  the  gradual  build-up  of 
circulation  about  the  airfoil,  and  should  be 
combined  with  the  momentum-integral  method  to 
produce  a  better  airfoil  pitching  model. 

b)  Equations  (54)  and  (55)  resulted  from  making 

fewer  assumptions  than  those  used  in  the  acutal 
computer  program.  A  procedure  on  how  to  use  them 
is  given  on  page  22;  the  computer  program  should  be 
updated  accordingly.  As  explained  previously,  this 
method  should  be  more  accurate:  not  only  were 

fewer  assumptions  made,  in  its  development,  but  the 
closure  equation  (5:13)  is  not  required  for  its 
solution.  It  could  also  be  used  for  the  study  of 
the  time  rate  of  change  of  the  boundary- layer 
thickness  due  to  pitching  motion. 


Further  work  needs  to  be  done  to  resolve  the  question  of 
whether  mass  ingestion  and  induced  camber  model  the  same 
physical  phenomenon.  If  they  do,  then  the  relaxation 
coefficient  used  in  Lawrence  mass  ingestion  model  is  not  four 
but  something  larger.  Perhaps  this  relaxation  coefficient 
should  be  a  function  of  position  along  the  airfoil  (i.e., 
larger  coefficients  could  be  rationalized  near  the  leading 
edge  where  relaxation  distances  are  smaller). 
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Appendix  A 

Derivation  of  the  Unsteady,  Non-Inertial  Euler  Equation 

For  a  non-inertial  reference  frame.  Shames  shows 
(7:141-144)  that: 

dF  -  dm  [R  +  2tox  Vx yZ  +  to  x  r  +  to  x(to  x  r)]  =  D_  /Jm  t?  \  (A.l) 


DtXyz  (<Jm  Vxyz^ 


For  the  case  where  pressure  P  is  the  only  surface  force 
(i.e.  no  viscous  force)  and  the  weight  of  the  fluid  in  the 
c.v.  is  the  only  body  force,  we  have 


dF  =  -g(VZ)  p  dvol  -  Vp  dvol 


(A. 2) 


If  we  substitute  (A. 2)  in  (A.l)  and  divide  thru  by 
pdvol  =  dm,  we  get: 


'I  VP  "  9Vz  “  [  R  +  2<o  x  Vxyz  +  w  X  r  +  wx(u)  X  r)] 


=  L 

Dt 


'xyz 


xyz 

=  (VXyz  •  v)  Vxyz  +  IV 


xyz 


xyz 
9txyz 


(A. 3) 


This  is  the  general  non- ine r t x a  1 ,  unsteady  Euler 
equation.  We  now  apply  (A. 3)  at  the  edge  of  the  boundary 
layer,  in  the  c.v.  shown  in  Figure  3  of  the  main  text. 

We  take  only  the  x-component  of  (A. 3)  and  neglect 
gravity  forces,  (note  here  that  even  if  we  had  not  neglected 
these  in  both  the  momentum  equation  and  Euler  equation,  they 
would  cancel  out)  we  get: 


-_1  __3P  -  {R  +  2co  XV  +  a)  x  r  +  u  x  (to  x  r)} 
p  3x  yZ 
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=  u3u  +  v3u  +  3u_ 

3x  3y  3t 


(A.4) 


A.  . 


From  the  main  text,  the  x-component  of  the  expression  in 
brackets  is,  for  our  specific  c.v.: 

R(asin(J>  +  a2cos<}>)  -  2av  -  ay  +  a2x 

And  therefore: 

-  _1  j3P  -  R(asin<f>  +  a2cos<}))  -  2av  -  ay  +  a2x  =  UoU+vau+3u  (A. 5) 

p  ax  ax  ay  at 

Applying  (A. 5)  at  the  edge  of  the  boundary  layer  and 
noting  that  8  Ue/3y  =  0 ,  we  get: 

-  _ 1.  3P_  _  ue  3 Up  +  3 Up  +  R(asin<}>  +  acos<|))  +  2ave  +  afi-a2x  (A. 6) 

p  dx  a  x  at 

(A. 6)  is  the  Euler  eguation  applied  at  the  edge  of  our 
particular  non-inertial ,  unsteady  control  volume.  As  3P/3x 
is  by  definition  a  function  of  x  only,  its  value  is  the  same 
across  the  height  of  the  boundary  layer  for  any  given  x.  The 
value  given  by  (A. 6)  can  therefore  be  substituted  in  the 
momentum  equation. 
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Appendix  B 


Unsteady  Separation 

For  a  steady  boundary  layer,  separation  occurs  when  the 
shear  stress  at  the  wall  is  zero;  combining  this  with  the  no¬ 
slip  condition,  we  see  that  the  separation  criteria  is: 

at  y  =  0  U  =  0 

3u/9  y  =  0 

For  a  moving  wall,  however,  the  no-slip  condition 
becomes 

at  y  =  0  U  =  UWALL 

In  this  case  the  Moore-Root-Sears  (MRS)  criterion  for 
separation  is  that: 

U  =  0  (B. la) 

9u/  3y  =  0  (B.  lb) 

at  some  point  in  the  interior  of  the  boundary  layer  (8:427), 
as  seen  by  an  observer  in  the  inertial  reference  frame.  In 
this  Appendix,  we  combined  the  RMS  criterion  with  the 
momentum-integral  method  to  find  the  value  of  A  at  separation 
for  our  pitching  airfoil. 

In  our  development  of  the  von  Karman-Poh lhausen  method, 
we  have  assumed  a  velocity  profile  of  the  form: 

~  =  A  +  Bn  +  Cn2  +  Dn3  +  Eg4  (B,2) 


53 


and  found  that 


A  =  0  (B. 3a) 

B  =  A  +  12  (B.3b) 

6 

C  =  -  J_  (B.3c) 

2 

D  =  A  -  4  (B.3d) 

2 

E  =  -  A  +  6  (B .  3e) 

6 


Substituting  (B.3)  into  (B.2),  and  applying  (B.lb),  we  find 
that  for  separation  to  occur,  the  velocity  profile  must 
satisfy : 


+  3ff  -  y 


~  V  +  (2  + 


(B .  4 ) 


where  A  s  is  the  value  of  the  velocity  profile  shape 
parameter  A  at  separation.  The  cubic  polynomial  (B.4)  can  be 


solved  analytically  (14:9).  By  applying  the  CRC  method  of 


solution,  we  find  that  (B.4)  has  three  real  roots: 


These  are  then  the  value  of  n  at  which  9  ( U /  )/  =  0  or 

equivalently  3u/3y  =  0.  The  double  root  occurs  atri  =  1,  i.e. 
at  the  edge  of  the  boundary  layer  where  u  =  ue  and  therefore 
3u/  3y  =  0  as  normal.  The  first  (B.5a)  must  then  be  the 


solution  for  the  moving  wall  separation.  We  note  that  it 
includes  the  steady  wall  solution  because  if  we  let  q  ^  be 
zero,  we  obtain  Ag  =  -12  which  is  the  steady  wall  separation 
criterion.  Substituting  (B.5a)  back  into  (B.2)  ,  we  get 

Uc  =  (As  +  12) 2  (9  A?  -  136Ac  +  528)  (B.6) 

Uoo  "  51 2 (As  -6) 3 

From  the  inertial  reference  U  =  0  and  again  we  recover 
A  =  -12  as  the  separation  criterion.  From  the  non-inertial 
reference  frame,  however,  Us  =  -  UWALL  or,  for  our  pitching 
airfoil  problem: 

n  •  n  .  ,  (B.7a) 

Us  =  -aRsin<{> 

Uoo  =  Ue  (B.7b) 

Combining  (B.6)  and  (B.7),  we  get: 

(As  +  12)2  (9A#  -  136AS  +  528)  =  -  512aRsin(}>  =  DK  (B.8) 

(As  -  6  )3  ue 

For  given  airfoil  and  test  condition,  (B.8)  can  be 
solved  for  As.  Figure  B.l  is  a  plot  of  the  solution  of 
equation  (B.8)  as  a  function  of  DK.  When  the  value  of  A  s  is 
reached  at  the  quarter  chord,  separation  is  assumed  to 
occur . 
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Repeats  program  for  next  test  conditions  from  FL0WIN2 
I F( ALPHA . LE. 900) GO  TO  125 

FORMAT <4X,F6.3,2(4X,F10.3),4X,F7.3,4X,F7.4,4X,F8.4,2(4X,E9.3)) 
FORMAT( ///////) 

FORMATC 1H1 , "BOUNDARY-LAYER  PARAMETERS  FOR  u , F6 . 2 , “FT/SEC" / ) 
FORMAT ( “  INITIAL  ANGLE  OF  ATTACK:  “,F6.3,“  DEGREES"/) 

FORMAT( / "  FINAL  ANGLE  OF  ATTACK:  “,F6.3,“  DEGREES”/) 

FORMATS  PITCH  RATE:  ",F7.3,M  DEGREES/SEC" / ) 

FORMATC 6X, "XOC“ , 10X, "U" , 1 1 X , “ DUDS" , 9X , "DUDT" , 8X, “FK" , 

9X, "RK" , 1 IX, “Z" , 1 IX, "DZDS*  /) 

FORMATC  PITCH  PARAMETER:  ",F7.5/> 

FORMAT < "  K  AT  THE  QUARTER-CHORD:  ”,F8.4/) 

FORMATC  STALLING  COEFF  <DK>:  ”,F8.4/) 

FORMATC "chord=  “ ,F8 . 4, 10X, "angle  for  0.25r  ”,F8.4////) 

FORMAT ( "  TIME  TO  REACH  THE  QUARTER-") 

FORMATC  CHORD  FROM  THE  STAGNATION  POINT:  .  ?7.5,“  SEC"/) 
FORMAT ( "  UWALL/UE  AT  THE  QUARTER-CHORD:  “,E9.3/> 

FORMAT< "a irf o i 1  geometry  amu=  " , f6 . 3, 3x, "bmu=  ",f6.3/) 

FORMAT( "rotat i on  point  in  fraction  of  chord  (+  after  0.5c", 

"-  before):  ",f9.3 /) 

STOP 

END 


SUBROUTINE  U ( ANGLE , RADIUS , CON , El , UI NF , AMU , BMU , ALPHA , UU ) 
COMPLEX  CMPLX, Z, El , DZETA , W 

Function  of  this  subroutine  is  to  to  compute  the  local 
value  of  velocity  on  a  Joukowski  airfoil  using  complex 
potential  flow  theory. 

ANGLA= ANGLE- ALPHA 
X=RADIUS*COS< ANGLA*CON) 

Y=RADIUS*SIN< ANGLA*CON) 

Z=CMPLX(X, Y) 

W=UINF*C  < 1 . , 0 . ) - ( RADI  US* *2 ) /Z**2+ 
C2.*EI*RADIUS*SINCALPHA*C0N) ) /Z) 

X=X+AMU 

Y=Y+BMU 

Z  changed  to  represent  values  of  coordinates  used  in 
the  transformation  equation. 

Z=CMPLX(X, Y) 

DZETA=  <Z**2- 1 . )/Z**2 
UU=CABS< W) /CABS(DZETA) 

RETURN 

END 


J-m 


■ 


f 


I 


I 
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CALL  U< ANGLE , RADI  US, CON, El ,UINF, AMU, BMU, ALPH 1 ,U2C> 

CALL  NIU( ANGLE, RADIUS , CON , AMU , BMU, ADOT2 ,  X 1 ,  Y  1 , XROT 
,U2C,U2,RAS2,RAC2,DK> 

Compute  the  unsteady  velocity  gradient. 

DUDT= ( U2-U 1 ) /DELT 
ANGLE = ANGLE- 0 . 0 1 
ANGLE  1= ANGLE- 0. 01 

CALL  U( ANGLE  1 , RADIUS, CON, El ,UINF, AMU, BMU, ALPH1 ,U2C) 

CALL  DS< ANGLE  1 , RADIUS , CON , AMU , BMU , X2 , Y2 ) 

CALL  NI U( ANGLE  1 , RADI  US, CON, AMU, BMU, ADOT2,X2,Y2, XROT 
,U2C,U2, RAS2 , RAC2 , DK ) 

ANGLE0=ANGLE+0 . 0 1 

CALL  U(ANGLE0, RADIUS, CON, El , UI NF , AMU , BMU , ALPH 1 ,U0C> 

CALL  DS ( ANGLEO , RADIUS , CON , AMU , BMU , XO , YO ) 

CALL  NIU (ANGLEO, RADIUS, CON , AMU , BMU , AD0T2 , XO, YO , XROT 
, UOC , UO, RASO , RACO, DK  > 

CALL  U (ANGLE, RADI  US, CON , El , UI NF , AMU , BMU , ALPH 1 , U 1 C) 

CALL  DS( ANGLE, RADIUS , CON , AMU , BMU , X 1 , Y 1 ) 

CALL  NIU(ANGLE, RADIUS, CON, AMU, BMU, AD0T2, XI ,Y1 , XROT 
, U 1C, U 1 , RAS 1 , RAC  1 , DK) 

Compute  arc  length  and  velocity  gradient. 

DS1 =(SQRT( (X1-X0)**2+(Y1-Y0)**2>> /CHORD 

DS2  =  (SQRT( (X2-X1 )**2  +  ( Y2-Y 1 >**2>  > /CHORD 

DSS=DS 1 +DS2 

DUDS= (U2-U0) /DSS 

UDUDS=U 1 *DUDS 

X0C=(X1+XLE) /CHORD 

Stop  the  computation  at  the  quarter-chord,  but  only  when 
it  is  reached  on  the  upper  surface. 

IF((XOC. GE. 0.250) .AND. (XOC . GT . OLDXOC) >  GO  TO  25 
OLDXOC=XOC 

IF( N . LT . 250)  GO  TO  20 
N=0 

WRITE( 16, 1 )X0C,U1 ,DUDS,DUDT,FK,RK,ZZ,DZDS 
CONTINUE 

WRITE( 16, 1 )X0C,U1 ,DUDS,DUDT,FK,RK,ZZ,DZDS 

WRI TE( 1 6 , 45 ) ALPH 1 

WRI TE( 16,55)PITCH 

WRITE( 16, 60 )RK 

WRITE( 16,80) 

WRITE( 16,81 )TIME 
UMIN=RAS1 /UI 
WRITE( 16,85)UMIN 


I 

s 


DS1=<SQRT< (XI -X0)AA2+(Y1-Y0) **2) ) /CHORD 
DS2=(SQRT(<X2-X1 )AA2+<Y2-Y1 )**2>> /CHORD 

* 

*  Compute  the  arc  length  and  the  velocity  gradient. 

* 

DSS=DS 1 +DS2 
DUDS= ( U2-U0) /DSS 
X0C=(X1+XLE) /CHORD 
IFCN.LT.50)  GO  TO  10 
N  =  0 

WRITE( 16, 1 )X0C,U1 ,DUDS,DUDT,FK,RK,ZZ,DZDS 
10  CONTINUE 

* 

AD0T=AD0T 1 
N=0 

RAS 1=0.0 
RAC  1  =  0.0 
DO  20  J=K 1 , K2 

A 

*  Function  of  this  loop  is  to  compute  the  behavior 

*  of  the  boundary  layer  as  it  is  subjected  to  a 

*  pitching  motion. 

A 

N=N+ 1 

* 

*  Compute  the  pertinent  boundary  layer  parameters. 

* 

ZZ=DZDS*DS 1 +ZZ 

RK=ZZ* ( DUDS*  < 1 .+RAS1/U1 )+DUDT/Ul) 

CALL  POHL(RK.RLAMDA) 

DEL2  =  37 . / 3 1 5 . -RLAMDA/945 . - C  RLAMDAa A2 ) /9072 . 

FK  =  2 . ADEL2 A (  2 . - . 3683ARLAMDA+ . 01 04aRLAMDAaa2+ 

+  (RLAMDA**3) /4536 . ) 

F1K=< . 3-RLAMDA/ 120.) /DEL 2 
RMI C=4 . 0 

* 

* 

*  RMIC  is  the  mass  injection  constant. 

* 

A 

DZDS= ( FK+ (4.+F1K) aZZaDUDT /U 1  +  C  4 . +2 . aFF 1 K- 2 . /DEL2  >  AZZ 
+  aRAS 1 aDUDS/U 1 -2 . ARMI CaZZaRAC 1 /DEL2  > /U 1 

A 

a  Compute  the  time  increment  for  a  particle  to 

a  travel  from  point  Ci)  to  point  (i+l>. 

A 

DELT=CH0RDADS1/U1 
TIME=TIME+DELT 
DALPHA=DELTAADOT 
ANGLE=ANGLE+DALPHA 
ALPH l =ALPH 1 +DALPHA 
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* 

DUDS=(U2-U0) /(DS1+DS2) 

* 

*  Second  derivative  of  velocity  computed  using  a 

*  Taylor's  series  expansion. 

* 

D2UDS2=<U2-2.*Ul+UO)/<<DSl+DS2>/2. >**2 

•k 

*  Enter  initial  boundary  layer  parameters. 

* 

RLAMDA=7 . 052 
RK=0 . 0770 
FK  =  0. 0 

DZDS=-0. 0652*D2UDS2/ (DUDS* *2 ) 

ZZ=RK/DUDS 

* 

N=50 

ANGLE=2*ALPHA+THETA-0 . 0 1 
XOC=(XO+XLE) /CHORD 

WRITE< 16, 1 )XOC,UO,DUDS,D2UDS2,FK,RK,ZZ,DZDS 

* 

AD0T=0. 0 
DO  10  J= 1 , K 

* 

*  Function  of  this  loop  is  to  compute  boundary  layer 

*  parameters  at  stagnation  point,  allowing  the 

*  boundary  layer  to  steady-out  before  subjecting  it 

*  to  a  pitching  motion. 

* 

N=N+ 1 

* 

*  Compute  pertinent  boundary  layer  parameters. 

* 

ZZ=DZDS*DS 1 +ZZ 
RK=ZZ*DUDS 
FK= . 47-6 . *RK 
DZDS=FK/U 1 

* 

DELT=CHORD*DS 1 /U 1 
TIME=TIME+DELT 

CALL  U  <  ANGLE , RAD I  US , CON , E I , U I NF , AMU , BMU , ALPHA , U2 ) 
DUDT= (U2-U 1 ) /DELT 
ANGLE= ANGLE-0. 01 
ANGLE 1 = ANGLE- 0 .01 

CALL  U< ANGLE  1 , RADIUS, CON, El, UINF, AMU, BMU, ALPHA, U2> 
CALL  DS  <  ANGLE  1 , RAD I  US , CON , AMU , BMU , X  2 , Y  2 ) 
ANGLE0=ANGLE+0 . 0 1 

CALL  U < ANGLEO , RADIUS , CON , El , U I NF , AMU , BMU , ALPHA , UO ) 
CALL  DS < ANGLEO , RADIUS , CON , AMU , BMU , XO , YO ) 

CALL  U < ANGLE , RADIUS , CON , El , U I NF , AMU , BMU , ALPHA , U 1 ) 
CALL  DS< ANGLE, RADI US, CON, AMU, BMU, XI ,  Yl> 
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XN= 1 80. 

XM= ( XP+XN )  /  2 . 

CALL  DS<XM , RADI  US , CON , AMU , BMU , XC , YC ) 

IF<ABS<XP-XN) .LT. 0.00001 )  GO  TO  11 
FFXM=  <  < XC+XLE) /CHORD) -0 . 250 
IF(FXM.LT.O.O)  GO  TO  8 
IF(FXM.EQ.O.O)  GO  TO  11 
IF(FXM.GT.O.O)  GO  TO  9 

XN  =  XM 
GO  TO  7 
XP  =  XM 
GO  TO  7 
CONTINUE 
ANGLQC=XM 

CALL  UtANGLQC, RADIUS, CON, El ,UINF, AMU, BMU, ALPHA, UC) 

CALL  NIU< ANGLQC, RADI  US , CON, AMU , BMU, AD0T2 , XC, YC , XROT , UC, 
UNI ,RAS2,RAC2,DK) 

DK=-5 1 2 . *AD0T2*DK 

PI TCH=AD0T1 *C0N*0. 5*CH0RD/UINF 


K=  100 
K 1 =K+ 1 

K2= ( ALPHA+ 1 50) *  1 00+K 

WRITE( 16, 12) 

WRI TE< 16, 30)UINF 
WRI TE< 16, 40) ALPHA 
WRITE< 1 6 , 50) ADOT 1 
WRI TE< 16, 89)AMU, BMU 
WRITE< 16, 88 )XROT 
WRITE( 16, 65)CH0RD, ANGLQC 
WRITE( 16,66)DK 
WRITE< 16,52) 

ANGLE=2*ALPHA+THETA 

CALL  U ( ANGLE , RAD I  US , CON , E I , UI NF , AMU , BMU , ALPHA , UO ) 

CALL  DS< ANGLE, RADIUS , CON, AMU , BMU , XO, YO) 

ANGLE= ANGLE-0. 01 

CALL  U  <  ANGLE , RADI  US , CON , E I , UI NF , AMU , BMU , ALPHA , U 1 ) 

CALL  DS< ANGLE, RADI  US , CON, AMU, BMU, X 1 , Y 1 ) 

ANGLE=ANGLE-0 . 0 1 

CALL  U ( ANGLE , RADIUS , CON , El , U I NF , AMU , BMU , ALPHA , U 2 ) 

CALL  DS( ANGLE, RADIUS, CON, AMU, BMU, X2,Y2) 

DS2  =  <SQRT<  <X2-X1)**2+(Y2-Y1)**2)) /CHORD 
DS 1  =  <SQRT( <X1-X0)**2  +  <Y1-Y0)**2)) /CHORD 

Stagnation  point  velocity  gradient  computed  using  a 
forward  difference  method;  all  other  velocity  gradients 
computed  usiing  central  difference  method. 


*  BETA  =  local  slope  of  airfoil  surface. 

*  RAS1,RAC1  =  functions  (of  geometry)  due  to  pitching  motion. 

*  RLAMDA,RK  =  Pohlhausen  shape  parameters. 

*  FK , ZZ, DZDS , DEL2 , F 1 K  =  functions  of  shape  parameters. 

*  XROT  =  Location  of  the  rotation  in  percent  chord,  with 

*  respect  to  the  mid-chord  position.  A  negative  is 

*  interpreted  as  forward  of  the  mid-chord  position. 

*  TH1  =  Angle  on  the  circle  mapping  to  the  training  edge 

*  TH2  =  Angle  on  the  circle  mapping  to  the  leading  edge 

*  DK  =  The  quarter-chord  separation  parameter  accounting 

*  for  the  MRS  separation  condition  of  a  moving  wall 

*  RMIC  =  The  mass  ingestion  constant 

* 

PROGRAM  POHLIO 
COMPLEX  El 

OPEN  < 15,FILE=’FLOWIN2' ) 

REWIND  15 

OPEN  < 16,FILE='FLOWOUT2' ) 

REWIND  16 
El  =  <  0 . ,  1  .  ) 

125  READ  ( 15,*,END=100)  ALPHA , ADOT 1 , UI NF, AMU , BMU , XROT 
RADI US=SQRT ( ( 1 . -AMU) **2+BMU**2> 

OLDXOC= 1 . 

ALPH 1 =ALPHA 
CON=3. 1415927/180. 

AD0T2=AD0T 1 *CON 

TH 1 =- <  AT AN ( BMU / < 1 . -AMU ) ) ) /CON 

TH2= 180. -TH 1 

THETA= 180. 

TIME=0. 0 

* 

*  Compute  chord  length 

* 

CALL  DS < TH 2 , RADIUS , CON , AMU , BMU , XLE , YLE ) 

CALL  DS< TH 1 , RADIUS, CON, AMU , BMU , XTE, YTE) 

C=SQRT( ( YLE-YTE)**2+(XLE+XTE>**2) 

CHORD=C 

4  TH2=TH2+0 . 02 

CALL  DS < TH 2 , RADIUS , CON , AMU , BMU , XC , YC ) 

C=SQRT ( ( YTE- YC ) **2+<XTE-XC)**2) 

IFCC.GT. CHORD)  THEN 
CHORD=C 
GO  TO  4 
Et  SE 
GO  TO  6 
END  IF 

6  ANGLE=TH2-0. 01 

CALL  DS< ANGLE, RADI  US, CON, AMU, BMU, XLE, YLE) 

XLE=ABS(XLE> 

* 

*  Find  the  angle  mapping  to  0.25c 
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Appendix  D 


Computer  Program  POHL  KP  84S 


POHL  KP  84S:  pitching  airfoil  with  mass  introduction. 

NOTE  mass  ingestion  may  be  turned  off  by  assingning 
0.0  to  the  mass  ingestion  constant  RMIC,  on  line  235 


ALLAIRE,  Andre  J.S. 
GAE-84S 


This  program  is  adapted  from  program  ‘‘P0HL2",  Docken, 
GAE-82D,  and  "P0HL6 ",  Lawrence ,  GAE-83D.  It  will  be  used 
to  continue  analysis  work  in  the  area  of  dynamic  stall 
effects.  The  representative  airfoil  is  a  15%  Joukowski 
pitching  at  a  constant  rate  about  its  mid-chord  in  a 
steady  freestream;  however,  any  thickness  or  camber 
ratio,  as  well  as  location  of  rotation  point,  may  be 
specified.  If  desired,  mass  may  be  injected  into  the 
boundary  layer  from  the  freestream  to  model  the  non- 
Newtonian  motion  of  the  airfoil. 


Symbology  used  within  this  program  is  as  follows: 


El  =  imaginary  unit;  i.e.,  square  root  of  -1. 

RADIUS  =  radius  of  circular  cylinder. 

AMU,BMU=of f set  distance  of  center  of  circular  cylinder. 
ALPHA  =  angle  of  attack,  in  degrees. 

AD0T1  =  pitch  rate,  in  degrees  per  second. 

PITCH  =  non-dimensional  pitch  rate. 

CHORD  =  chord-length  of  Joukowski  airfoil. 

ANGLE  =  radial  of  a  given  point  on  the  circular  cylinder. 
UINF  =  freestream  velocity,  in  feet  per  second. 

CON  =  conversion  factor,  degrees  to  radians. 

X,Y,U  =  coordinates  and  potential  flow  velocity  on  airfoil. 
Z,W  =  complex  coordinates  and  velocity. 

UMIN  =  ratio  of  wall  velocity  to  potential  flow  velocity. 
TIME  =  cumulative  time. 

DELT  =  increment  of  time. 

K,N  =  integer  counters. 

DSS  =  increment  of  distance  on  the  airfoil  surface. 

DUDS  =  spatial  partial  derivative  of  velocity. 

DUDT  =  partial  derivative  of  velocity  wrt  time. 

XOC  =  chord  position. 
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to  consider  separately  the  upper  and  lower  surface.  From 
Figure  C4,  we  see  that  for  the  upper  surface: 


a  =  Si  +  63  (C.  12) 

r  =  /  f2(x)  +  (x  -  xrP'  (C.13) 

tan  61  =  f 1  (x)  (C.  14) 

tan  82  =  f (x)  (C •  15 ) 

(x  -  xr  ) 

84  =  H  +82-33  (C.16) 

2 

and  from  the  law  of  sines: 

^ _  =  tur _  (C .  17) 

sin  84  sin  83 


for  six  equations  in  six  unknowns.  Unfortunately,  the  system 
is  transcendanta 1  and  cannot  be  solved  analytically. 
However,  one  can  see  from  the  figure  C4,  that  the  head  of  the 
vector  uir  is  shifted  to  the  left  as  a  result  of  thickness, 
and  we  get  a  higher  local  angle  of  attack  than  for  the  flat 
plate,  which  results  therefore  in  the  creation  of  more 
induced  camber.  This  additional  camber  could  explain  the 
higher  AC£  predicted  by  Tupper  for  the  J015  airfoil. 
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Jumper  (15)  has  shown  that  for  the  case  of  a  parabolic 


camber  distribution,  the  zero  lift  angle  of  attack  is 


a£o  =  ‘2ZC/c  (C,8) 

where  Zc/c  is  the  camber  ratio.  Using  this,  our  pitching 
plate  has: 

alo  =  -«ND/2  (C>9) 

For  a  thin  flat  plate  at  small  angle  of  attack,  the 
theoretical  lift  coefficient  is  (16:2-8) 

C  £  =  Zna  ( C .  1 0 ) 

The  rotation  has  therefore  induced  an  increase  in  lift 
coefficient  of: 

AC£  =  Ha^jQ  (C .  11) 

In  his  trailing  edge  vortex  method,  Tupper  found, 
empirically/'  that  for  a  flat  plate  A  =  3.14  a  ND  our 

analytical  solution  confirms  his  result  with  remarkable 
accuracy.  Furthermore,  for  a  Joukowki  J015  airfoil,  Tupper 
found  (6:42)  a  jump  in  c^  of  3.47  ocND.  Without  solving  the 
solution  analytically,  the  next  section  provides  an  argument 
showing  thickness  may  increase  the  value  of  the  induced 
camber . 

Thick  Airfoil 


The  same  procedure  can  be  repeated  for  the  case  of  a 
general  thick  and  cambered  airfoil.  In  this  case,  we  need 


To  find  the  equation  of  the  equivalent  cambered  plate,  we 
simply  need  to  solve  the  differential  equation  (C.3).  As  w, 
xr  and  U  are  not  a  function  of  x,  the  solution  of  equation 
(C.3)  is: 

f(x)  =  -  0)  x2  +  x.cox+k 
2U°°  U°° 

The  integration  constant  k,  is  evaluated  by  imposing  f (x)  =  0 
at  x  =  xLE  =  0  (leading  edge).  Then  k  =  0  and 

f(x)  =  -  u)  /x2  -  2xrx\  (C .  4) 

2U°°  V  / 

If  we  non-dimensionalize  with  respect  to  the  chord: 


f(x/c)  =  -  aND  x2  "  2x  xr 
c2  ~TZ 

where 

aND  =  uiC 
2Ua> 


(C .  5 ) 


This  equation  shows  that  the  induced  camber  is  parabolic.  By 
equating  the  first  derivative  of  (C.5)  to  zero,  we  find  that 
the  maximum  of  the  camber  function  occurs  for  (x/c) 
(xr/c),  i.e.  at  the  rotation  point.  There, 


f (xr/C)  =  a(xr/C)2 

Finally,  we  note  that 


(C .  6 ) 


f(l)  =  2(xr/C)  -  1 


(C .  7) 


For  the  special  case  where  xr  =  0.5c,  f  ( 1 )  =  0  and  the 
induced  camber  ratio  is  a  ND/ 4 . 
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Appendix  C 
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Camber  Induced  by  Rotation 

When  an  airfoil  rotates  in  a  steady  flow,  the  local 
angle  of  attack  changes  along  the  surface  of  the  airfoil,  as 
can  be  inferred  from  the  pitching  flat  plate  of  figure  Cl. 
It  is  then  possible  to  model  a  pitching  airfoil  by  a  non¬ 
pitching  one  having  some  camber  induced  by  the  pitching 
motion.  The  amount  of  camber  equivalent  to  a  given  pitching 
rate  is  computed  in  this  Appendix.  We  consider  first  the 
simple  case  of  a  flat  plate  instantaneously  at  zero  angle  of 
attack  before  generalizing  to  a  thick  airfoil. 

Pitching  Flat  Plate 

From  figure  C2,  we  see  that  the  local  angle  of  attack 
seen  by  a  pitching  flat  plate  ,  for  a  general  point  x  on 
the  plate, is  given  by: 

tan  a(x)  =  to(x  -  xr)  (C.l) 

In  the  case  of  a  non-pitching  but  cambered  plate, 

if  f(x)  represents  the  locus  of  points  of  the  camber  in  the 
same  flow,  its  local  angle  of  attack  is  given  by: 

tan  a(x)  =  -f ’ (x)  (C. 2) 

Equating  (C.l)  and  (C.2)  gives: 

f 1  (x)  =  -  o)(x  -  xr)  (C.3) 
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* 

SUBROUTINE  DS( ANGLE, RADI  US , CON , AMU , BMU ,  X ,  Y  > 

COMPLEX  CMPLX , 2 

* 

*  Function  of  this  subroutine  is  to  compute  the  pos- 

*  iton  of  a  point  on  the  airfoil,  given  its  position 

*  on  the  mapping  circle. 

* 

X=RADI US*COS( ANGLE*CON  >  +AMU 
Y=RADIUS*SIN( ANGLE*CON ) +BMU 
Z=CMPLX(X, Y) 

Z=Z+1 . /Z 
X=REAL(Z> 

Y=AIMAG(Z) 

RETURN 

END 

* 

* 

SUBROUTINE  POHL< RK , RLAMDA ) 

* 

*  Function  of  this  subroutine  is  to  compute  the  value 

*  of  the  separation  parameter,  lamda,  giiven  a  value 

*  of  K,  as  computed  in  the  main  program. 

* 

RK 1 =- . 160 
RK2=- . 112 
RK3=0 . 00 
RK4=0 . 06 
RK5=0. 076 
RK6=0 . 086 
RK7=0. 0949 

IF(RK.LE.RKl)  GO  TO  10 
IF(RK.LE.RK2)  GO  TO  20 
IF( RK . LE . RK3)  GO  TO  30 
IFCRK.LE.RK4)  GO  TO  40 
IF(RK.LE.RK5)  GO  TO  50 
IFCRK.LE.RK6)  GO  TO  60 
I FC RK . GT . RK7 )  GO  TO  70 

* 

RLAMDA = .0149**2-(RK-0.08)**2 
RLAMDA= 12.-100. *SQRT<  RLAMDA ) 

RETURN 

* 

10  RLAMDA =  <  2 . / . 0 1 2 ) *RK+ 14.0 
RETURN 

* 

20  RLAMDA =  <  4 . / . 044 ) *RK+2 .  1  8 
RETURN 

* 

30  RLAMDA=  < 10- /. 14) *RK 
RETURN 
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* 

40  RLAMDA=83. 33*RK 

RETURN 

* 

50  RLAMDA=- 1 .9+1 1 5 . *RK 

RETURN 

* 

60  RLAMDA=-6. 54+176. *RK 

RETURN 

* 

70  RLAMDA= 1 2 . 

RETURN 

END 

* 

* 

SUBROUTINE  N I U ( ANGLEC , RADIUS , CON , AMU , BMU , AD0T2 , XC , YC , 

+  XROT, UC, UNI , RASF, RACF, DK> 

* 

*  Function  of  this  subroutine  is  to  compute  the  value  of  U  in 

*  the  non-inertial  reference  frame,  and  the  values  of  the 

*  required  functions  of  geometry,  as  a  function  of  the  airfoil 

*  shape  and  the  location  of  the  rotation  point.  The  value  of 

*  the  quarter-chord  separation  parameter  DK  is  also  computed 

*  this  subroutine. 

* 

R=SQRT(XC**2+YC**2 ) 

BET 1 =ATAN (YC/<-XC>) 

ANGLEM=ANGLEC+ . 1 
ANGLEP=ANGLEC- . 1 

CALL  DS < ANGLEM , RADIUS , CON , AMU , BMU , XM , YM > 

CALL  DS < ANGLEP , RADIUS , CON , AMU , BMU , XP , YP ) 

DY=ABS( YP- YM) 

DX=ABS(XP-XM) 

Q= 1 000 . *DX 
IFCDY.GT.Q)THEN 
BETA=1 .5707963 
ELSE 

BETA=BET1+ATAN(DY/DX> 

END  IF 

RP=SQRT<  YC**2  +  (XR0T-XC >  **2 ) 

ANG=ABS<  <RP**2+R**2-XROT**2  > / <  2 . *R*RP) ) 

IF(ANG.GE. 1 . 0  )  THEN 

DELTA  1=0.0 

ELSE 

DELTA  1 =AC0S<  ANG) 

END  IF 
R=RP 

I F  <  XROT . GE . 0 . 0 ) THEN 
BETA=BETA-DELTA 1 
ELSE 

BETA=BETA+DELTA1 
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